function [f, df] = PSFest(k,u_hat,b_hat,u_u,u_b,R)
%PSF_EST Uses minFunc to estimate a PSF kernel, k
% The function outputs the value of the objective function, f and a vector
%  of the objective derivative, df.
% Author: Rose Rustowicz, rose.rustowicz@gmail.com
% Date: 16 March 2018

% Set constants
% lambda = 0; 
mu = 10;
gamma = 0.1;
eta = (2*R+1)*(2*R+1);
k_size = sqrt(size(k(:),1));

k = reshape(k, [k_size, k_size]);
dx = [0,  0, 0;
    0, -1, 1;
    0, 0, 0];
dy = [0,  0, 0;
    0, -1, 0;
    0,  1, 0];



% --------- Compute f, the Expected value of k. ---------- %
% First term
f1 = mu .* sum(sum(k(:).^2));
df1 = 2 * mu .* k(:);

% Second term
f2 = norm(u_hat * k(:) - b_hat).^2;
% df2 = 2 *(u_hat' * u_hat * k(:)) - (2*u_hat' * b_hat);
df2 = 2 *(u_u * k(:)) - (2*u_b);

% Third term
kx = conv2(k, dx,'same');
ky = conv2(k, dy,'same');
f3 = gamma * eta  * ((sum(sum(kx(:).^2))) + sum(sum(ky(:).^2)));
dxt = circshift(ifft2( conj(fft2(circshift(dx, [-1 -1]))) ),[1 1]);
dyt = circshift(ifft2( conj(fft2(circshift(dy, [-1 -1]))) ),[1 1]);
dxk = conv2(dxt,dx);
dyk = conv2(dyt,dy);
fft2_k = fft2(k);
dfx = 2 * gamma * eta * ifft2(fft2_k .* psf2otf(dxk, size(k)));
dfy = 2 * gamma * eta  * ifft2(fft2_k .* psf2otf(dyk, size(k)));
df3 = dfx + dfy;
df3 = df3(:);

f = f1 + f2 + f3;
df = df1 + df2 + df3;

end

